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ABSTRACT 

A binary supermassive black hole loses energy via ejection of stars in a galactic nucleus, until 
emission of gravitational waves becomes strong enough to induce rapid coalescence. Evolution via the 
gravitational slingshot requires that stars be continuously supplied to the binary, and it is known that 
in spherical galaxies the reservoir of such stars is quickly depleted, leading to stalling of the binary 
at parsec-scale separations. Recent A^-body simulations of galaxy mergers and isolated nonspherical 
galaxies suggest that this stalling may not occur in less idealized systems. However, it remains unclear 
to what degree these conclusions are affected by collisional relaxation, which is much stronger in the 
numerical simulations than in real galaxies. In this study, we present a novel Monte Carlo method that 
can efficiently deal with both collisional and collisionless dynamics, and with galaxy models having 
arbitrary shapes. We show that without relaxation, the final-parsec problem may be overcome only 
in triaxial galaxies. Axisymmetry is not enough, but even a moderate departure from axisymmetry is 
sufficient to keep the binary shrinking. We find that the binary hardening rate is always substantially 
lower than the maximum possible, “full-loss-cone” rate, and that it decreases with time, but that 
stellar-dynamical interactions are nevertheless able to drive the binary to coalescence on a timescale 
< 1 Gyr in any triaxial galaxy. 

Subject headings: galaxies: evolution — galaxies: kinematics and dynamics — galaxies: nuclei 


1. INTRODUCTION 

Binary supermassive black holes (SBHs) are naturally 
formed in galaxy mergers, if both merging galaxies con¬ 
tain a central SBH. The heavy objects quickly sink to the 
center of the merger remnant due to dynamical friction 
and form a binary system. Subsequent evolution of the 
binary is driven by interaction with stars in the galactic 
nucleus, which are ejected by the slingshot mechanism 
(|Saslaw et al.lIlQTll if they arrive within a distance < a 
from the binary center of mass, where a is the semi¬ 
major axis of the binary orbit. As a result, the orbit 
shrinks (hardens), and if the binary becomes sufficiently 
hard, emission of gravitational waves (GWs) becomes 
the main source of energy loss, rapidly bringing the two 
SBHs to coalescence. The efficacy of this process, how¬ 
ever, depends crucially on the supply of stars into the loss 
cone (the low-angular-momentum region of phase space 
in which stellar orbits can approach the binary), and if 
the reservoir is depleted, the binary nearly stalls, or at 
least its shrinking tim escale may become m uch longer 
than the Hubble time (|Begelman et al.lll980f) . 

In idealized spherical galaxies, the only guaranteed 
mechanism of loss-cone repopulation is two-body relax¬ 
ation. The problem of feeding stars into the loss cone 
of the binary has much in common with the similar 
problem for a single SBH, which was extensively stud¬ 
ied i n the 1970s in the c ontext of spherical system s 
(e.g. iFrank fc ReeiJ 119761 iLightman fc Shapirol il977l l. 
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These papers identified an important distinction between 
empty- and full-loss-cone regimes: in the former, the flux 
of stars is inversely proportional to the relaxation time 
and depends on the size of the loss cone only logarith¬ 
mically, while in the latter the relaxation is so efficient 
that the supply of stars into the loss cone becomes inde¬ 
pendent of the relaxation time and proportional to the 
size of the loss cone. Relaxation times in real galaxies 
are so long that they are nearly always in the empty- 
loss-cone regime. The conjectured stalling of the binary 
evolution has been labeled the “final-parsec problem” 
(jMilosavlievic fc Merrittll2003a|) . 

On the other hand, in non-spherical galactic potentials 
the angular momentum of stars changes not only by two- 
body relaxation, but also by large-scale torques (jMerrittl 
120131 Ghapter 4). The orbital structure of non-spherical 
galaxies is rather diverse, and in triaxial galaxies there 
exist an entire class of centrophilic orbits, (e.g. box or 
pyramid orbits) that may attain arbitrarily low values of 
angular momentum without any relaxation. These orbits 
were identified as a promising mechanism of loss-cone 


repopulation (e.g. INorman & Silklll983 

; IMerritt & PoonI 

120041: iHollev-Bockelmann & Sigurdsson 

I2OO6II. Similarly, 


in an axisymmetric potential only one component of an¬ 
gular momentum is conserved, and the number of stars 
that can enter the loss cone is much larger than in the 
spherical case (iMagorrian fc Tremaind 119991 120021 : 

iVasiliev fc MerrittI 120131) . although not as large as in a 
triaxial system. 

Numerical simulations of the evolution of binary SBHs 
face an important difficulty: since the number of parti- 
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cles in a typical simulation (< 10®) is several orders 
of magnitude smaller than the number of stars AA in a 
real galaxy, it is necessary to properly understand the 
scaling laws. The rate of collisional evolution is inversely 
proportional to the relaxation time, which scales roughly 
as N/logN, while collisionless effects are essentially in¬ 
dependent of N. Even if only the collisional effects play 
a role in the dynamics, the hardening rate scales dif¬ 
ferently in the empty- and full-loss-cone regimes. For 
small N the system is in the latter regime and the hard¬ 
ening rate is nearly independent of N, while for large 
N the hardening rate should drop with N. Early stud¬ 
ies were restricted to rather low values of N and hence 
did not find any A^-dependence of the hardening rate 


(e.g. Quinlan fc HernauisdllQQTtlMilosavlievic fc MerrittI 


1200111 . while in more recent simulations of isolated 
spherical systems with larger N t he hardening rate 
was f ound to decline w i th N (e.g. [Makino fc Funatol 
12004 iBerczik et al.ll200^ iMerritt et al.l 1200711 . although 
less steeply than the A^“^ dependence expected in the 
empty-loss-cone regime. On the oth er hand, simulations 
that considered isolat ed triaxial (IBerczik et al.l 120061: 
Berentzen et al.l 1200911 or axisymmetric (iKhan et al.l 


201311 systems, or started from a merger of two galax¬ 


ies, which need not result in a spherical model, typically 
find no dependence of hardening rate on N. This has 
been interpreted as a sign that the binary remains in the 
full-loss-cone regime due to efficient reshuffling of angular 
momenta of orbits by non-spherical torque s. 

In a previous paper (|Vasiliev et al.l[2014L hereafter Pa¬ 
per I), we reconsidered binary hardening in isolated 
galaxies with different geometries (spherical, axisymmet¬ 
ric and triaxial) using high-resolution A'^-body simula¬ 
tions with N up to 10®. Somewhat surprisingly, we 
found that the hardening rates do depend on N in all 
three cases, although they decline less rapidly with N 
for non-spherical models. With the exception of the 
highest-A^ integrations, there was almost no difference 
between axisymmetric and triaxial models. We also ex¬ 
plored the possible contribution to the hardening rate 
from collisionless effects, by analyzing the properties of 
orbits in our models and estimating the draining rates 
of centrophilic orbits. We concluded that with presently 
accessible values of N it is difficult to disentangle colli¬ 
sional and collisionless effects in loss-cone repopulation. 
In real galaxies, however, collisional effects are expected 
to play a much smaller role, so that it is hard to draw firm 
conclusions about the evolution of binary SBHs in non- 
spherical galaxies based on conventional A^-body simula¬ 
tions. 

In this study, we return to the topic and present a 
novel Monte Carlo method that can be used to model the 
evolution of galaxies hosting binary SBHs. Our new al¬ 
gorithm contains an adjustable rate of relaxation, which 
can even be set to zero, yielding the collisionless limit. 
We demonstrate that in this limit, triaxial systems have 
enough centrophilic orbits to maintain an adequate sup¬ 
ply of stars into the loss cone, although the binary hard¬ 
ening rate is not as high as in the full-loss-cone regime 
and slowly declines with time. Nevertheless, for all rea¬ 
sonable values of the parameters, the coalescence time 
is shorter than the Hubble time. We therefore conclude 
that the final-parsec problem is solved by triaxiality even 


in a purely collisionless stellar system. By contrast, in 
axisymmetric and spherical galaxies the hardening rate 
rapidly drops in the absence of relaxation, meaning that 
in most galaxies the binary would never merge. Crucially, 
collisional relaxation in conventional N^-body simulations 
overwhelms the depletion of the loss cone and completely 
changes the long-term behavior of the binary. 

In section we describe the Monte Carlo method 
used in this study and compare it to previous similar 
approaches, while in section [3] we validate the method 
against a large suite of conventional A^-body simulations 
with N < 2 X 10®. Section [4] is devoted to a detailed 
study of the evolution of isolated galaxy models, con¬ 
structed initially as equilibrium configurations in spher¬ 
ical, axisymmetric and triaxial geometries. We present 
the results of Monte Carlo simulations and illustrate the 
trends found in long-term evolution with simple analyt¬ 
ical arguments. Based on these results, we compute co¬ 
alescence times for binary SBHs as a function of galaxy 
structure and initial parameters of the binary orbit. In 
section[5]we conduct N-body simulations of mergers and 
compare them to Monte Carlo models. Finally, in sec¬ 
tion [6l we summarize our important results and compare 
them wi th previo us work on the final-parsec problem. 

As in iPaoer D . we focus on purely stellar-dynamical 
processes. For a broader picture, includi ng the ef f ects o f 
gas-dynamical torques, see Se c tion 8 .4 of lMerrit^ (120131 1 
or the recent review of iColoil (l2014| l. Some preliminary 
results f r om the work presented here were described in 
IVasilievI (l2014b[ l. 


2. METHOD 
2.1. Definitions 

We consider the evolution of a binary SBH composed 
of two point masses, mi and m 2 , which are on a Keple- 
rian orbit with semimajor axis a and eccentricity e. The 
total mass of the binary, Mbin = mi -|- m 2 , is a small 
fraction (10“^ — 10“®) of the total mass of surround¬ 
ing stellar distribution. The mass ratio of the binary is 
q = m 2 /mi < 1. 

A star passing at a distance ;< a from the binary un¬ 
dergoes a complex scattering interaction and is ejected 
with a typical velocity ~ \/ GM^^/a (the characteristic 
orbital velocity of the binary); for a hard binary this is 
higher than the average velocity of the stellar population, 
thus the star carries away energy and angular momentum 
from the binary. The precise definition of a hard binary 
varies among different studies; here we adopt that a bi¬ 
nary is hard if its semimajor axis is smaller than 


Oh = 


4(1+ g)- 




( 1 ) 


where, in turn, the radius of influence rinfl is defined as 
the radius enclosing the mass in stars equal to twice the 
total mass of the binary. This definition depends on the 
evolutionary phase, as the slingshot process reduces the 
density of stars in the galactic nucleus. The most rapid 
depletion occurs just after the binary formation, and af¬ 
ter the binary becomes hard the depletion slows down 
considerably. For consistency with merger simulations, 
in which it is not possible to assign any particular value 
to rinfl before the two galactic nuclei have merged, we 
adopted the convention to evaluate the influence radius 
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after the hard binary has formed (see iMerritt fc Szelll 
120061 for a discussion on different definitions of rinfl). 

Another important quantity is the hardening rate S = 
d(l/a)/dt. In what follows, we will frequently compare 
S with the reference value S'f^n that would occur if the 
distribution of stars in phase space were not affected by 
the presence of the binary - in other words, if the loss 
cone was “full”. Unfortunately this value also does not 
have a commonly accepted definition. For instance, if the 
stars were uniformly distributed in space, with density p 
and with isotropically directed velocities all of the same 
magnitude v, then the hardening rate can be expressed 
as 


5*uniform — ^1 5 

V 


( 2 ) 


where the dimensionless coefficient Hi can be mea¬ 
sured from scatte ring experiments (e.g. iQuinlanI 119961 : 
ISesana et al.ll2006ll . and has a value Hi ~ 18 in the hard¬ 
binary limit. A more widely used definition applies to 
systems with a uniform density and a Maxwellian veloc¬ 
ity distribution with dispersion cr: 


f{x,v) = pf{v) = 


(27ro-2)3/2 


exp - 


2ct2 


.Gp 


»^Maxwell — / dv AtTV /(u) (S'uniform ("c) — H , (3) 

Jo cr 


where H = Hiy/2/Tr 14.5. 

The assumption of a uniform-density background is 
clearly an oversimplification, and a more robustly de¬ 
fined quantity is obtained by averaging the single-velocity 
hardening rate over the actual distribution function 
of stars f{x,v) in the entire galaxy. If we assume that 
the latter is isotropic (i.e., depends only on the energy 
U, but not on the angular momentum), then 




2 TJ 

dv 4?™^ f{E) —-— = A-k HiG 


dEf{E). 

( 4 ) 


Here 4>o is the depth of potential well of the stellar 
cusp (excluding the potential of the SBH), so that the 
integration includes energies corresponding to stars that 
are unbound to the binary, but still bound to the entire 
galaxy (this is a rather ad hoc convention, but it is justi¬ 
fied by the fact that few stars remain bound to the binary 
for any significant time after it becomes hard). Another 
derivation of this quantity fr om a som ewhat different per¬ 
spective can be fou nd e.g. in Paper J . Equations (7)-(9), 
or in Section 8.3 of IMerrittI ( 201311 . Note that the above 
value has no relation to the total mass of stars in the 
galaxy, which is given by / dE f{E) g{E), where g{E) is 
the density of states and rapidly rises with E. 

Yet another way to define a reference hardening rate 
is to substitute the density and velocity dispersion eva l- 
uated at nnS into Equation ([3]) (jSesana fc Khanir2015l l. 
This definition has the advantage of being easily com¬ 
puted from quantities that are accessible to an A^-body 
snapshot, without the need to calculate the distribution 
function; on the other hand, it reflects only the local 
parameters near the galaxy center and not the entire 
population of stars that can interact with the binary. 


Since the velocity dispersion can itself be expressed as 
cr2 ^ GMbin/nnfl, we may write 

5infl = A(GMbi„)i/2rr^/", (5) 

which agrees quantitatively with the if Here the dimen¬ 
sionless factor A can be chosen to match the previous 
definition (|1]) ; its value depends weakly on galaxy struc¬ 
ture and lies in the range 3-5. 

In what follows, we will refer to either Siso or S'inS 
as the “full-loss-cone hardening rate”. We stress that 
there is no fundamental reason to expect that the actual 
hardening rate will be close to ^fun, since the distribution 
of stars at low angular momenta is not isotropic, and even 
the value of Sfun need not remain constant in time - it 
merely serves as a useful reference value. 

The loss of energy and the evolution of eccentricity 
due to the emission of G Ws are given by the following 
expressions (IPeterslll964H : 


d{l/a) 

dt 

fie) 

de^ 

dt 


1/a 


Tqw = 




(1 + 9)= 


Tqw 

(l-e2)+2 


64 G3M3 


fie), (6) 


l + ie2. 


17 

96*^ 


( 7 ) 


q 26^(304-b 121e2) 
A’a'^ (1 + q/^ 15(1 — e2)3/2 


Throughout the paper, we will mostly present the re¬ 
sults in dimensionless A7-body units, with the mass and 
the scale radius of the galaxy model both set to unity. 
The models may be scaled to a given galaxy using any 
two out of three fundamental scales (length, time and 
mass). To simplify the discussion, we reduce this free¬ 
dom to a one-parameter family of galaxies in which the 
length and mass scales are related through the — a 
relation in the following form (e.g. IMerritt et al.l I2009L 
Fig.l2): 

HnS = T’o X (M,/10® Mq)'^, ro = 30 pc, k = 0.56. 

( 9 ) 


2.2. Monte Carlo code 

The Monte Carlo metho d used in this work is an ex¬ 
tension of the Raga code (lVasilievll2015ll . We follow the 
motion of N particles in the combined potential of the 
stellar distribution, $*(r), and the two point masses or¬ 
biting each other, centered at origin. The orbit of the 
binary is assumed to be Keplerian, aligned in the x — y 
plane, and elongated in x direction; we make no attempt 
to follow either the change of the orbital plane, which we 
know to be small from the A7-body simulations (although 
it could be quit e significant ove r long timescale in triax- 
ial galaxies, see lCui fc Yull201^ . or the precession of its 
periapsis, which we assume is not particularly important 
for dynamics. The evolution is broken down into many 
small intervals of time (episodes) of duration Tepi Tbin; 
during each episode we keep the stellar potential and the 
parameters of the binary orbit unchanged. 

Each particle is moving independently from the oth¬ 
ers under the influence of three forces: the gradient 
of the smooth static stellar potential, represented as 
a spherical-harmonic expansion with spline-interpolated 
coefficients as functions of radius; the time-dependent 
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force from the binary; and random velocity kicks that 
model the effect of two-body encounters for a system 
composed of a certain number of stars, which does 
not need to be related to the number of particles in the 
simulation, and can even be set to infinity (thus switch¬ 
ing off the two-body relaxation). The velocity perturba¬ 
tions are computed from the local (position-dependent) 
drift and diffusion coefficients, using the standard for- 
malism from the relaxation theory (e.g. iMerrittI I2013L 
Chapter 5) and an isotropic approximation for the distri¬ 
bution function of stars (i.e. with dependence on energy 
only). Thus our method descends from the Spitzer’s for- 
mulation of Monte Carlo approach (|SDitzer fc HartlfigTlf) 
and does not rely on orbit-averaging, as other con tempo¬ 
rary M onte Carlo codes that use the formalism of lHenonl 

Cml). 

During the encounter of a particle with the binary, de¬ 
fined as the time when the distance of the particle from 
origin is less than 10a, we record the changes in energy 
and angular momentum of the particle that arise due 
to the motion in time-dependent potential of the binary 
(not including the perturbations from relaxation or the 
torque from the non-spherical stellar potential). After 
each particle’s trajectory has been calculated over the 
entire episode, we sum up these changes for all particles 
that experienced one or more encounters with the binary, 
and change the binary’s energy and angular momentum 
by the same amount with the opposite sign, while keeping 
its orbital phase unchanged from the end of the episode 
to the start of the next one. Thus the stars do not exert 
force on the black holes directly, but we use conservation 
laws to impose the changes to the binary orbit param¬ 
eters. The binary orbit parameters are optionally mod¬ 
ified after each episode acco rding to the expressions for 
gravitational-wave emission (I6ia . The stellar potential 
and the diffusion coefficients that account for two-body 
relaxation are also updated at the end of each episode, 
reflecting the changes in the stellar distribution in the 
course of evolution (here the most important effect is the 
gradual erosion of the central stellar cusp). 

2.3. Comparison with similar approaches 

It is instructive to compare our method with the 
other schemes used for an approximate treatment of the 
joint evolution of the binary and the stellar distribution. 
iQuinlan fc HernanistI (119971 1 developed a program SCF- 
BDY which combine s elements from the self-co nsistent 
field (SCF) method (|Hernauist fc Ostrike]Hll99^ with a 
direct integration of black hole-star interactio ns. The 
two b lack holes were integrated using NB0DY2 (lAarsethl 
[19^ - a direct-summation code with neighbor scheme, 
adaptive timestep, and optional two-body regularization, 
taking into account the forces from all stars individually. 
The motion of star particles, on the other hand, was 
computed in the gravitational potential of the two black 
holes and the smooth potential o f the other stars, rep¬ 
resen ted by a basis-set expansion (iHernauist k. Ostrikei] 
Il992fl . Their approach is quite similar to ours, with the 
following differences: (a) the binary’s center of mass is 
not fixed at origin, although the center of stellar poten¬ 
tial expansion is; (b) the star particles exert “real” forces 
on the black holes, which lead directly to the changes 
in the binary orbit, instead of relying on the Newton’s 
third law as in our method; (c) two-body relaxation 


is not modeled. They used a rather small update in¬ 
terval for the stellar potential, which could lead to ar¬ 
tificial relax ation due to fluctuations i n the expansion 
coefhcieiits (IHernauist fc BarnesI 119901: iWeinbera 119961 
lSellwoodll2015h : in our code we recompute the potential 
less frequently and use Agamp ^ 1 sampling points from 
each particle’s trajectory stored during each episode, 
to further reduce discreteness noise. Nevertheless, any 
finite-A systern is not entirely free of relaxation (e.g. 
iWeinber^ 119981 1. so our simulations without explicitly 
added two-body relaxation should be regarded as an up¬ 
per limit for the evolut ion rate expected in truly coll ision- 
less systems. Finally, IQuinlan fc HernauistI (|1997D only 
dealt with nearly spherical systems, although in princi¬ 
ple their scheme can equally well work for non-spherical 
systems after switching on the corresponding terms in 
t he spherical-harmo n ic exp ansion (as in our method). 

iHemsendorf et aH (120021 1 used a similar approach in 
their code EuroStar: the two black holes and a sub¬ 
set of star particles with angular momenta lower than a 
certain threshold are integrated with a direct-summation 
method, adapted from Aarseth’s NBODy6 code, while the 
rest of the star particles are followed in the smooth field 
b y a modified yersioii of the SCF method. 

iSesana et al.l (|2006l I2008D used a large suite of three- 
body scattering experiments to derive the expressions for 
the rate of ha rdening and ecc entricity growth, extending 
the results of iQuinliml (1199611 to a wider rm ge of mass 
ratios and eccentricities. ISesana et all (j2007[l and lSesaii^ 
(j2010ll developed a hybrid approach for the evolution of 
the binary, combining the analytical fits to these scatter¬ 
ing experiments with a time-dependent model for the loss 
cone draining and repopulation. They did not simulate 
the effect of two-body relaxation and Brownian motion 
of the binary center of mass explicitly, but included some 
prescriptions to take it into account. The changes in the 
st ellar potential d u e to e jection of stars were ignored. 

iMeiron fc Laori (|2012ll introduced another scheme 
based on the conservation laws to find the changes in the 
binary’s orbit. The motion of stars in the static spher¬ 
ically symmetric potential of the stellar cusp plus the 
time-dependent potential of the binary is followed over a 
short interval of time, after which the changes in the total 
energy and angular momentum of all stars are translated 
into the forces that should be applied to the black holes. 
Unlike our method, in which we assume a Keplerian orbit 
for the binary and adjust its parameters on a timescale 
much longer than the orbital period, their approach re¬ 
quires rather short update intervals to follow the motion 
of the black holes directly. 

Compared to the previo us studies, our approach most 
clos ely resembles that of IQuinlan fc HernanistI ()1997r i 
and IHemsendorf et aH (1200211 . in that we also use the 
spherical harmonic expansion technique to represent the 
smooth potential of the stellar cluster, which is itself a 
nearly collisionless system, while considering the inter¬ 
action between stars and the binary as a succession of 
three-body scattering events, as in the series of papers by 
Sesana et al. The treatment of the binary is somewhat 
more approximate in our method - we assume its center- 
of-mass to be at rest, and ignore the changes in its orbital 
plane and orientation of the line of nodes. The evolution 
of its most important parameters - binding energy and 
eccentricity - is deduced from the orbits of stars using 
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conser vation laws, in a similar manner to lMeiron fc La^ 
(|2012ll . but using much longer update intervals, spanning 
many orbital periods. The most important new feature 
of our method with respect to the previous studies is 
the ability to efficiently model the two-body relaxation 
with adjustable magnitude, allowing us both to compare 
our Monte Carlo simulations with direct iV-body simu¬ 
lations (using the same relaxation rate), and to extend 
them into nearly collisionless regime, by switching off the 
relaxation. 

Fokker-Planck methods have also been applied 
to this problem, always in the spherical geometry 
(iMilosavlievic fc Merrittl [20031)1 : iMerritt fc Wand 120051 : 
iMerritt et al.ll2007ll . The Fokker-Planck approach has 
the advantage of a potentially much finer resolution 
of the loss-cone region in {E, L) space. A disadvan¬ 
tage is the need to orbit-average the diffusion coefh- 
cients. Our Monte Carlo algorithm avoids that limi¬ 
tation; on the other hand, the orbit-averaged approx¬ 
imation only breaks down near the loss cone, and at 
least in the spherical geometry, a boundary-layer treat¬ 
ment of the loss cone is available that is based on the 
local (not orbit-averaged) Fokker-Planck equation and 
which automatically acc ounts for empty vs. full loss cones 
(|Cohn fc Kulsrudlll978ll . The secondary slingshot is also 
more naturally accounted for in a Monte Carlo algorithm, 
although it has also been implem ented in Fokker-Planck 
treatments (|Merritt et al.l l2007fl . The primary advan¬ 
tage of our method over Fokker-Planck algorithms is the 
ability to model stellar systems regardless of their shape, 
including the chaotic orbits that arise in non-spherical 
geometries. Limitations of our approach include neglect 
of Brownian motion (justified to some extent below in 
Section 14.31) and the assumption of near-equilibrium of 
the stellar distribution, which means that it is not ex¬ 
pected to perform well in highly dynamic situations such 
as mergers. It is also worth noting that the method scales 
linearly with N and already at N = 10^ outperforms the 
GPU-accelerated direct-summation code by a factor of 
20 . 

3. TESTS 

We first checked that the results of our calculations 
do not depend on the technical parameters such as the 
number of radial grid points in the spherical-harmonic 
expansion, the update interval, or the type of orbit inte¬ 
grator and its accuracy parameters, provided that they 
are set to reasonable values. 

We also verified that the conservation-law method pro- 
du ces the s ame results as the scat t ering experiments 
of iQuinl^ (|1996ll and iSesana et al.l (|2006ll which used 
a small but non-zero mass of the incoming star, and 
computed its effect on the binary orbit directly. In 
particular, the hardening rates and eccentricity growth 
ra tes were fouiid to b e comparable to Figures 3 and 4 
of ISesana et al.l (|2006ll for a uniform-density isothermal 
background population of stars. 

Next we compared the Monte Carlo code with a large 
suite of conventional A^-body integrations of isolated 
galaxy models with various shapes and initial parame¬ 
ters for the binary. These were similar to the simulations 
we used in iPaoer but covered a wider range of eccen¬ 
tricities, mass ratios, and densit y profiles . The initial 
density profile of stars follows the iDehnenI (Il993f) model 


with the inner cusp slope 7=1 (our default model) and 
7 = 2, in three different geometries: spherical, oblate 
axisymmetric (axis ratio 1 : 0.8) and triaxial (axis ratios 
I : 0.9 : 0.8 - our default model, - and 1 : 0.8 : 0.6). 
These are constructed to be in equilibrium with a cen¬ 
tral SBH of mass M, = O.OI of the total stellar mass, and 
have a nearly isotropic velocity distribution. The mass 
M, is split between two components of the binary, and 
the two SBHs are initially placed at a separation 0.2 (for 
7 = I) or 0.02 (for 7 = 2), slightly larger than the radius 
of influence Q We considered two values for the mass ra¬ 
tio - q = 1 (equal-mass binary) and q = 1/9. Using dH), 
our 7 = l,q = 1 models can be scaled to real galaxies 
so that one length unit equals 150 pc x (M,/10® 
and one time unit is 0.27 Myr x (M,/10® 

The A^-body simulations were p erformed with th e 
direct-summation code c^iGRAPEch (|Harfst et al.ll200^ . 
using the SAPPORO library for GPU acceleration 
(|Gaburov et al.l[200^ . This code employs chain regular¬ 
ization to accurately follow the motion of the binary and 
the stars interacting with it, in exact Newtonian gravity 
(no softening); we used a very small softening e = 10“® 
for particles outside the chain. Thanks to the use of the 
chain, the relative error in energy was typically below 
10“® at the later stages of evolution, when both massive 
particles are included in the chain. 

The initial conditions for Monte Garlo models were 
taken from A^-body simulations at the time when the 
binary has just become hard, reaching the semimajor 
axis flh. The relaxation rate is a free parameter in the 
Monte Garlo method, defined by the number of stars in 
the target system AA and the Coulomb logarithm In A. 
We have measured the changes in energies and angular 
momenta of particles in the spherical A^-body simula¬ 
tions, and compared them t o the expected diffusion rate 
as a function of energy Isee lVasilievll20I5l . figure 1); the 
value of A = 0.02AA provided the best agreement for 
10® < Ni, < 10®. The prefactor in the Coulomb loga¬ 
rithm is the only free parameter in the code that could be 
assigned at will (apart from technical parameters such as 
the number of terms in the potential expansion), and we 
have adopted the above definition to match the diffusion 
rate, not the hardening rate of the binary or anything 
else; thus if the other aspects of evolution agree between 
Monte Carlo and A^-body simulations, this demonstrates 
the predictive power of the Monte Carlo approach. 

We compared the Monte Carlo and A^-body simula¬ 
tions using a number of criteria. The most important 
are the binary parameters - semimajor axis and eccen¬ 
tricity. It should be noted that individual simulations 
have a considerable scatter in the hardening rates and 
the evolution of eccentricity, especially at low N. Never¬ 
theless, statistically the agreement between Monte Carlo 
and A-body simulations is very good for a wide range 
of parameters (A, initial eccentricity, geometry, binary 
mass ratio), see Figure [T] for a few examples. We have 
also checked that the long-term behavior of Monte Carlo 
simulations with the same relaxation rate (set by A*) 

^ Recall that we measure ri^fj just after the hard binary has 
formed; for instance, in a 7 = 1 Dehnen model with a single SBH 
rinfl = 0.165, but after the binary becomes hard, the density cups 
is eroded and r^^fi increases to ~ 0.2 for the models with g = 1 , 
changing only slightly in the subsequent evolution. 
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Fig. 1.— Evolution of inverse semimajor axis as a function of 
time, for a few models with various N, initial eccentricity, binary 
mass ratio g, and geometry. Thinner lines are from A^-body simu¬ 
lations, and thicker ones are from Monte Carlo simulations. 



a^/a 



Fig. 2.— Evolution of binary eccentricity as a function of sep¬ 
aration, in Ai-body (top panel) and Monte Carlo (bottom panel) 
simulations. Shown are tracks for isolated systems with equal-mass 
binaries (solid red) and g = 1/9 (dashed green), and merger simu¬ 
lations (dot-dashed purple); the adopte d an alytic prescription for 
the evolution of eccentricity (Equation IIOII is shown by the blue 
lines. 

but with different number of particles N is similar; this 
allows to extrapolate our method to large /V* while still 
using a reasonable {N < 10®) number of particles. 

The dependence of hardening rate on eccentricity was 
found to be weak; if anything, models with high e evolved 
a little faster, in agreement with the results of scat¬ 


tering experiments of iSesana et al.l (|2006D . [Merritt et'aTI 
(|2007l) . but typically the difference was comparable to the 
scatter between individual runs. The evolution of eccen¬ 
tricity itself is in fact more important, since the energy 
losses due to GW emission depend strongly on e. Scat¬ 
tering experiments typically suggest a slow but steady 
growth of eccentricity, but in practice its evolution is 
more erratic, because individual interactions may both 
increase and reduce e (thus the outcome results from a 
slight imbalance between them) . Stars with smalle r an¬ 
gular momenta tend to reduce e ([Sesana e t al.||200 8D. and 
so do stars that corotate with the binary ([Iwasawa et all 
1201 It ISesana et aT]l 201 lD . In our simulations, the frac¬ 
tion of corotating and counterrotating stars is approx¬ 
imately equal, so the latter factor does not come into 
play, but the distribution of stars in angular momentum 
does depend on the details of loss-cone repopulation, so 
the former effect is quite important. The eccentricity it¬ 
self stayed roughly constant if it was initially small, and 
growed slightly if started from e > 0.5. Models with ini¬ 
tial e > 0.8 typically increased it to 0.9 < e < 0.95 in the 
course of evolution, especially in the unequal-mass case 
{q — 1/9). Nevertheless, neither in V-body nor in Monte 
Carlo simulations did we observe a rapid growth of eccen- 
trici ty to values > 0.99, found in some previo us studies 
(e.g. llwasawa et al.l20lHlMeiron &: Laoii2012D . although 
we did not consider systems with such high mass ratio 
as in the former study. 

The eccentricity growth is traditionally described with 
a dimensionless parameter K = de/dhi{l/a). Previ¬ 
ous studies have found that in the hard-binary limit, K 
reaches a maximum value of ^ 0.1 — 0.2 at e ~ 0.7, and 
drops to zero at e = 0 or 1. We adopt the following 
param etrization for K , whi ch is comparable to the find¬ 
ings of lQuinlanI (jl99(il ) and ISesana et al.l (| 200 fi[) : 


K = 


de 


Ae{l-ey 1 + 


c?ln(l/a) 
b = 0.6, oq = 0.2ah, A = 0.3 . 


ao 


-1 


( 10 ) 


Figure [2] compares the evolution of eccentricity based 
on the above equation (dotted lines) to the results of N- 
body and Monte Carlo simulations. Even though there 
is considerable scatter between runs, the overall trend is 
reasonably well described by the fitting formula (TTOl) . 

4. EVOLUTION OF ISOLATED GALAXY MODELS 


4.1. Long-term evolution of stellar-dynamical hardening 

rate 

We now consider the long-term evolution of binary 
SBH in galaxies with a realistically large number of stars 
/V*. Figure [3] shows the hardening rates S computed on 
the interval 50 < t < 100 for a series of 7 = 1 mod¬ 
els with different V*, shape and binary mass ratio. It 
appears that in both spherical and axisymmetric cases 
the hardening rates continue to drop with increasing iV*., 
but triaxial models tend to a nonzero limiting value of 
S as iV* —>• 00 . This was already suggested in iPaner II 
based on V-body simulations, but the number of parti¬ 
cles N < 10® was not enough to establish it clearly; ad¬ 
ditional simulations with N = 2 x 10® and Monte Carlo 
models support this conclusion. 

Figure [4] shows simulations extended to a much longer 
interval of time, approximately 0.4 Gyr for a model scaled 
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Fig. 3.— Hardening rates as functions of for different geome¬ 
tries and mass ratios, computed on the interval 50 < t < 100 from 
A^-body simulations (filled symbols, solid lines) and Monte Carlo 
simulations (open symbols, dashed lines). The agreement is quite 
good, but Monte Carlo models underestimate the hardening rate 
for Nir > 10 ® in the spherical case, presumably due to the neglect 
of Brownian motion. 

to a galaxy with Mbin = 10® Mq. For each geometry we 
show the simulation without relaxation (i.e. in the col¬ 
lisionless limit), as well as for a moderately small relax¬ 
ation rate (smaller than is achievable in 7V-body simu¬ 
lations). It is immediately clear that in the collisionless 
limit there is a striking difference between three geome¬ 
tries: in the spherical case, the binary stalls at a semima¬ 
jor axis barely smaller than Oh, and in the axisymmetric 
case it shrinks roughly a factor of ten below Ch, but ul¬ 
timately also stalls. By contrast, in the triaxial case the 
binary continues to shrink, although the hardening rate 
decreases with time. In the simulations that include re¬ 
laxation, however, the situation is quite different - the 
binary never stalls in any geometry. For the triaxial case, 
the hardening rate in the TV* = 5 x 10® simulation is 
already not much higher than in the collisionless limit, 
but the axisymmetric systems differ dramatically from 
their collisionless counterpart: in both TV* = 5 x 10® and 
N^, = 10® systems the hardening rate decreases at early 
times, as in the collisionless limit, but unlike the latter, 
it then attains a non-zero lower limit. 

We also considered a model with a stronger initial de¬ 
gree of triaxiality (axis ratio 1 : 0.8 : 0.6); initially it 
had a higher hardening rate, which then dropped to 
a level comparable to that of the less flattened model. 
While a higher hardening rate is naturally explained by 
a larger reservoir of chaotic orbits in a more triaxial 
model (see next section), its subsequent decrease is a 
result of the loss of triaxiality in the central parts of 
the model (Figure [S]). Interestingly, the model with a 
milder initial flattening retained its shape better. This 
evolution toward axisymmetry clearly occurs due to 
the binary: if we replace it with a single SBH, the 
shape stays nearly constant on much longer intervals 
of time than considered here. It is commonly assumed 
that even single SBHs ine vit ably destroy triaxial it y (e.g . 
iGerhard &: Binned 1198511 . iMerritt fc Quinlan~l i|1998 1 
simulated growth of black holes in triaxial models formed 



12 5 10 20 50 

aja 

Fig. 4.— Long-term evolution of binary hardness for the equal- 
mass case and three different geometries (Spherical, Axisymmetric 
and Triaxial). Shown are curves corresponding to different relax¬ 
ation rates (dash-dotted lines: A^* = 2 x 10®, 5 x 10® and TV* = 10®) 
and to the simulations with no relaxation (solid lines), dashed line 
is for the triaxial model with no relaxation and stronger initial 
flattening {y/x = 0.S,z/x = 0.6, while the other models have 
zjx = 0.8 and, in the triaxial case, yjx = 0.9). It is apparent 
that in the spherical and axisymmetric cases the binary separa¬ 
tion approaches an upper limit without relaxation, while in the 
triaxial case it continues to shrink. Also notable is that even a 
modest amount of relaxation keeps the binary from stalling, al¬ 
though the evolution rate could be quite low for a realistic number 
of stars (even for a very low-mass binary, TWbin = lO®TVff 0 , the 
model scaled to a real galaxy would have TV* = 10®). 

Top panel shows the time evolution of semimajor axis, and bot¬ 
tom panel shows the dependence of hardening rate on Uh/u. The 
full-loss-cone rate 0 is around 20 in our models, higher than any 
value measured in the simulations. Models with relaxation eventu¬ 
ally attain a nearly constant hardening rate, depending on TV* and 
geometry. In collisionless simulations, by contrast, the hardening 
rates keep decreasing, very mildly in the triaxial case and steeply 
in the axisymmetric case. Dotted lines show the as ymptot ic ex¬ 
pressions for the hardening rate in scale-free galaxies (|17I20|I . 

via cold collapse and found that central SBHs containing 
more than ^ 3% the mass of their host galaxies induced 
evolution toward axisymmetry. However, later studies 
have shown that this does not necessary happen in iso¬ 
lated stationary systems: using a more f l exible , orbit- 
superposition algorithm, iPoon fc MerrittI (I2002D found 
that self-consistent and apparently stable models could 
be constructed in which the SBH mass was a substantial 
fraction of the nuclear mass, a result confirmed in other 
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Fig. 5.— Evolution of shape of collisionless triaxial models. 
Dotted, dashed and dash-dotted lines show the axis ratio (y/x - 
top curves, blue and cyan, zjx — bottom curves, red and gold), 
measured at radii enclosing 10, 50 and 90% of total mass for the 
models with initial x:y\z = 1:0.9:0.8 (thinner lines) and 1:0.8:0.6 
(thicker lines with longer dashes). The latter model evolves much 
quicker toward axisymmetry in the central parts. 


studi ^es (iHollev-Bockelman et al.l 120021 : iPoon fc MerrittI 
12004) . iVasilievI 1 20151) demonstrated. using the same 
Monte Carlo code, that relaxation is the main driving 
force behind shape evolution: in the collisionless regime 
the shape remains nearly constant. The likely reason is 
that the diffusion of chaotic orbits, which tends to erase 
triaxiality, is greatly facilitated by the no ise from two- 
body relaxation (e.g. iKandrup et al.l[2000ll . The shape 
evolution seen in the collisionless simulations with a bi¬ 
nary SBH - as opposed to a single SBB - may be caused 
by resonant perturbations of c haotic orbits by the t ime- 
dependent gravitational field (|Kandrup et al.l 1200311 . al¬ 
though more work is needed to explore this effect. 

In addition, we performed simulations of models with 
a steep cusp (7 = 2 ) and otherwise the same parame¬ 
ters. They demonstrated a similar behavior, with only 
the triaxial model continuing to shrink indefinitely in the 
collisionless case, although the hardening rate was drop¬ 
ping with time more rapidly than in the 7=1 model. 
The loss of triaxiality was mild, much like the 7 = 1 
model with the same axis ratios (1:0.9:0.8). 

Before we proceed to a theoretical explanation of the 
hardening rates, we discuss the properties of orbits in our 
models. 


4.2. Properties of orbits interacting with the binary 

We have examined the types of orbits that bring stars 
into interaction with the binary. For each particle that 
arrives at a distance less than 10 a from the center of mass 
of the binary, we record its phase-space coordinates at its 
first approach. If a particle is scattered more than once 
by the binary, we only consider the first interaction, be¬ 
cause it is difficult to define unambiguously when one 
interaction ends and the next one begins. Then we com¬ 
pute a trajectory starting from these initial conditions 
in a smooth stellar potential plus a single point mass 
Mbin at the origin, for a time corresponding to 100 or¬ 
bital periods. This allows us to use powerful analysis 
tools applicable to orbits in smooth stationary poten¬ 
tials. To distinguish regular and chaotic orbits, we use 


the Lyapunov exponent, and to determine if an orbit is 
centrophi lic, we fo llow the algorithm described in the ap¬ 
pendix of lPaoerl . 

After the binary became sufficiently hard (a < 0.3 Oh), 
most of the orbits that interact with it are unbound to 
the binary (have energies higher than the depth of stel¬ 
lar potential well $ 0 ). In the non-spherical cases, almost 
all such orbits (> 90%) are found to be chaotic, and for 
the triaxial system, a similarly large fraction of them are 
truly centrophilic, i.e. may attain arbitrarily low values of 
angular momentum. There are no genuinely centrophilic 
orbits in axisymmetric systems, because the variation of 
angular momentum is bounded from below by its con¬ 
served z-component, but for chaotic orbits the average 
value of L is typically much larger than its minimum 
achievable value L^, thus they constitute a reservoir of 
“usable” orbits (loss region) that is much larger than the 
loss cone itself. 

Regarding the overall orbital structure of the model, 
the chaotic orbits are a minority (^ 10 %, depending on 
the degree of flattening and triaxiality). As discussed 
above, the shape of triaxial models gradually evolves to¬ 
ward axisymmetry from the inside out, but globally it 
remains sufficiently triaxial to support a substantial pop¬ 
ulation of centrophilic orbits - their fraction is roughly 
proportional to the deviation from axisymmetry, and 
their total mass is thus considerably larger than the mass 
of the binary. 

4.3. Theoretical models for the hardening rate evolution 

The evolution of binary in the collisionless limit can 
be qualitatively described with a rather simple model for 
the draining of the population of orbits that ca n interact 
with the binary - similar to the one presented in lPaDer~4 . 
but without a detailed analysis of properties of orbits in 
a particular simulation. 

We begin with the triaxial case, and assume that at 
each energy there is initially a fixed fraction 77 of chaotic 
centrophilic orbits, which occupy the low-angular- 
momentum region of phase space < riLl-j.^{E). To 
account for their gradual depletion, we introduce the 
fraction of surviving orbits f{E,t), so that the mass of 
chaotic orbits is 

Mch= / dEATT^T,,A{E)pLl^^{E)aE,t)f{E), (11) 
J 

where we have neglected the L-dependence of radial or¬ 
bital period Tj-ad- The hardening rate due to interaction 
with these low-angular-momentum orbits is given by a 
generalization of Equation |4j 

S{t) = = 47r iLi G [° dE f{E) f{E, t). (12) 

at 

On the other hand, the same interactions eject stars 
and decrease the fraction of surviving orbits. Equating 
the energy carried away by stars with the change in the 
binary’s binding energy gives 

== ( i) (13) 

2 a 2 \a J 

with = mim 2 /Mbin is the binary reduced mass (e.g. 
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lMerrittir2004 equation 25). Thus 


(iMch 

dt 


1 miTO 2 d 

3 /i ^ dt 



--MbinaS. 


(14) 


If we assume that the stars interacting with the bi¬ 
nary become unbound to the galaxy (i.e., neglect the 
secondary slingshot), then the draining of the loss region 
occurs independently at each energy, and we can write 
the following equation for the evolution of surviving frac¬ 
tion of stars: 


daE,t) ^ Hi 2GMbinQ(t) 

dt ’ ’ 6tt 'qLl.^^{E)T,^diE)' 


We take y = Oh/a as the new independent variable 
instead of t. Substituting (USD into the above equa¬ 
tion, we integrate it over y with the initial condition 
2 / = 1 , ^{E, y = 1) = 1 and obtain 


i{E,v) = exp 


HiGMun 




dy' 

y' s{y'y 


(16) 


Here h{y) is another unknown function related to S{y). 
Substituting the above expression back into (HU) , we 
can perform the integration over E, keeping h as a 
hxed parameter, and obtain an equation of the kind 
S{y) = g{h{y)) with some function g. Then we express 
h from the resulting equation as a function of S and y, 
and hnally differentiate h by y to obtain an ordinary dif¬ 
ferential equation for S{y). 

To illustrate this approach and acquire qualitative in¬ 
sight, we consider the asymptotic behavior of the system 
at large t for the case of a power-law density profile, 
p{r) oc r~^. Pure power-law profiles are unphysical in 
the sense that the total mass is infinite and the gravita¬ 
tional potential can take arbitrarily large positive values, 
but the contribution of stars at large radii to the hard¬ 
ening rate is negligible, so that all quantities of interest 
are finite. We choose to define the gravitational poten¬ 
tial of the stellar cusp so that its value at the center is 
$0 = 0 ; since we only consider stars that are not bound 
to the binary, their energies lie in the range 0 < E < oo. 
Moreover, we may neglect the potential of the SBHs, be¬ 
cause at late times most of the surviving loss-region stars 
are located at far larger distances than rinfl. Then all 
dynamical quantities have power-law dependence on en¬ 
ergy: f{E) oc £;-( 6 - 7 )/( 4 - 27 )^ L^irciE) OC i;( 4 - 7 )/( 4 - 27 ) ^ 
TradiE) OC ifT'/(4-27)^ Carrying out the integration in 
(fT^ . we obtain S{y) oc and ulti¬ 

mately get the asymptotic dependence of the hardening 
rate on a: 

2 + t 

S{a) oc [77/ln(ah/a)] . (17) 

In other words, the hardening rate due to depletion of 
centrophilic orbits drops with time (or a“^) very slowly, 
and has a moderate dependence on the fraction of chaotic 
orbits in the model rj. This expression is also valid in the 
limiting case of a singular isothermal sphere (7 = 2 ). 

We now follow a similar argument for the axisym- 
metric case. Here we again assume that orbits with 
< r]Ll-j.yE) are chaotic, but only those with < 


TABLE 1 

Long-term evolution of hardening rate in triaxial models. 

First column is the model type, second and third are the 
full-loss-cone hardening rate and the radius of hard binary, 
in model units; last two columns are the dimensionless pa¬ 
rameters of Equation 1211 computed as best-fit values to the 
hardening rate found in collisionless Monte Carlo simulations. 


model 

‘S'infi 

Oh 


V 

default (7 = 1 , q = 1 , 
x:y:z = 1 : 0 .9:0. 8 ) 

20 

0.0125 

0.38 

0.32 

<7 = 1/9 

30 

0.004 

0.33 

0.32 

x:y:z= 1 :0.8: 0.6 

20 

0.0125 

0.96 

0.51 

7 = 2 

1700 

0.0022 

0.25 

0.62 

merger (7 = 1), N=128k 

15 

0.017 

0.21 

0 

merger, N=256k 

15 

0.017 

0.43 

0.18 

merger, N=512k 

15 

0.017 

0.45 

0.25 

merger, N=1024k 

15 

0.017 

0.87 

0.47 


Llc = can interact with the binary. The 

mass of such “useful” chaotic orbits in the case of a hard 
enough binary, when Llc < is given by 

.0 

Mch= / dE T,.^yE)2^Lcirc{E)LLc{a) ^{E,t)f{E), 

J 

(18) 

and the expression for their depletion rate, analogous to 

(fT51) . is 


d^{E,t) _ Hi y2GMbina(t) 

dt ’ Utt 2^L,i,yE)T,ad{E)' ^ ’ 

The solution is obtained along the same lines as for 
the triaxial case, with the different definition of h(y) = 
f dy/{y/y S{y)). We again consider the asymptotical evo¬ 
lution described by dni), dUD for an idealized scale-free 
model and obtain 


S{a) oc 


{r]a/ah)^~^^ , 7 < 2, 

exp [-C'ah/(??a)] ,7 = 2 


( 20 ) 


Unlike the triaxial case, the hardening rate drops 
quickly with decreasing a. Moreover, in realistic non¬ 
scale-free systems the mass of chaotic orbits m that are 
both not yet depleted {^{E,t) ^ 1) and able to interact 
with the binary {L^ < Llc(o-)) is finite, and after some 
time drops below the mass of the binary itself, which 
means that the further evolution virtually ceases. 

These findings are beautifully illustrated by the results 
of collisionless Monte Carlo models in Figure IH bottom 
panel: the hardening rate in the triaxial case drops very 
gently with Oh/o and closely follows its asymptotic ex¬ 
pression dnD, while in the axisymmetric case it decays 
much faster, following (I20L and then drops nearly to zero 
when all chaotic orbits are depleted. Since our galaxy 
models are not scale-free and gradually lose the triaxial- 
ity, the hardening rate in the collisionless triaxial model 
is better described by a somewhat steeper dependence on 
a, namely S (x a'^ with v ranging from ^1/3 for mildly 
triaxial 7 = 1 models to > 1/2 for more strongly triaxial 
models or for steeper cusps (7 = 2 ), see TableHJ The fact 
that steeper cusps result in a stronger slowdown of hard¬ 
ening rate can be explained by a more rapidly declining 
dependence of the distribution function on energy: once 
the centrophilic orbits close to the binary are depleted, 
there are less available stars at larger distances in mod¬ 
els with steeper density profiles. Figure [ 6 ] illustrates the 
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Fig. 6 . — Illustration of the loss region and hardening rate evolution in various geometries (Spherical, Axisymmetric and Triaxial). 

Left panel shows the slice of the phase space (L, Lz) at a fixed energy E. The loss cone of the binary is the region L < Llc = \/2GM\^{^a\ 
stars in that region will interact and be ejected by the binary in one dynamical time. In non-spherical geometries, the extended loss region 
consists of stars with L < y/rjLdrc (more correctly, only a fraction of such stars are on chaotic orbits, but this does not qualitatively change 
the picture). In a triaxial geometry, these stars can wander anywhere in this region (shaded in blue and green) due to collisionless effects 
(non-spherical torques), and can eventually get into the loss cone proper (shaded in red). In the axisymmetric case they are restricted to 
lines of constant Lz and can only wander in vertical direction, thus the loss region is a strip Lz < Llc? ^ < \/V^ciTc (shaded in green). The 
population of stars inside the loss region is gradually depleted, and as the binary shrinks, Llc decreases. Dashed lines show the boundary 
of the loss cone and the axisymmetric loss region at a later moment of time; the volume of the loss region in the axisymmetric case also 
shrinks with the binary, and the number of stars in this region drops both due to its decreasing volume and decreasing fraction of surviving 
stars. 

Right panel illustrates the depletion of the loss region population as a function of energy and time. Bottom plot shows the fraction of 
surviving stars ^{E) at various moments of time. Initially ^ = 1 and it depletes faster at high binding energies; thinner curves correspond 
to later times. Top panel shows ^{E) multiplied by the distribution function; the integral under this curve gives the hardening rate at any 
given time, and initially is equal to the full-loss-cone hardening rate (Siso)- Thus it is natural that the hardening rate is always smaller 
than S'iso and decreases with time. 


above arguments about the size of the loss region and its 
depletion. 

Finally, we consider the effects of relaxation on the 
long-term behavior of the hardening rate. From the 
above discussion, it is clear that it may only matter for 
spherical and axisymmetric systems, because in the tri¬ 
axial case the draining is the main mechanism that keeps 
the loss cone filled. In general, the relaxation rate scales 
as N~^, but it does not trivially translate into the hard¬ 
ening rate because of several complications. First of all, 
the steady-state flux of stars into the loss cone has a dif¬ 
ferent dependence on the size of the loss cone Llc and the 
relaxation rate in the limits of empty and full loss cone 
regimes: in the former case, at a fixed energy it scales as 
[V* ln(Lcirc/I^Lc)]“^ and in the latter - as (iLc/^circ)^- 
Integrated over all energies with an appropriate bound¬ 
ary condition at each energy, the hardening rate has a 
weaker than N~^ scaling, which furthermore depends on 
the evolutionary stage: as a gets smaller, a larger fraction 
of total flux comes from the full-loss-cone region. Sec¬ 
ond, the flux of stars into the loss cone is actually higher 
than the steady-state models would predict, because ini¬ 
tially there are large gradients in the phase space: stars 
with L < Llc are absent while at larger L their dis¬ 
tribution is nearly unchanged. Time-dependent models 
for the loss cone repopulation (iMilosavlievic fc Merritt! 
I2003bf) predict a flux that is several times higher than 
the steady-state value at early times. 

Overall, we find that in the spherical case the hard¬ 
ening rate drops with V* rather mildly at fV* < 2 x 
10 ®, but for larger fV* approaches the asymptotic N~^ 
scaling. In iV-body simulations, however, it does not 


drop quite as fast at large iV*. One possible rea¬ 
son might be the wandering (Brownian motion) of 
the b inary (lOuinlan fc HernauisdllQQTl : iChatteriee et al.l 
l2003f) , which effectively increases the size of the loss cone 
from Llc = 2GMbina to = 2GMbinr'wand- The wan- 

derin g radius scales as r^and oc (m*/Mbin)®^^ oc 
(e.g. iMerrittl 120011) . thus for realistically large iV* it 
should remain below a for the most part of evolution, 
and will not substantially increase the hardening rate. 

Relaxation in the axisymmetric case occurs at the same 
rate as in the spherical system, but the effective size 
of the loss cone corresponds to the angular momentum 
of the chaotic region of the phase space y^Tcirc, which 
is then drained into the loss cone proper by the non- 
spherical torques and not by relaxation. This suggests 
that in the empty-loss-cone regime the steady-state flux 
is mo derately (a factor of few~) l arger t han in the spherical 
case (|Magorrian fc Tremainel 119991 : iVasiliev fc Merrill 
IMl), because it depends only logarithmically on the 
effective size of the loss region. Time-dependent flux is 
again higher than the stead y-state value, by larger fac tors 
than in spherical systems (jVasiliev fc Merritt! I201.H fig¬ 
ure 13, left panel). Giyen these uncertain complications, 
it is hard to derive more quantitative theoretical esti¬ 
mates for axisymmetric systems. The results of Monte 
Carlo simulations suggest that the hardening rate drops 

at least as V* in the range 10® < V* < 10®, and 
probably even steeper at larger iV*, which means that it 
is too slow for realistic galaxies to bring the binary to 
GW-dominated regime in a reasonable time. 
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4.4. Estimates of the coalescence time 

Motivated by the above arguments, we write the 
stellar-dynamical hardening rate S'* in an evolving model 
as 


S*(a) = AtSinfl (a/ah)'", (21) 

where the dimensionless coefficient /r < 1 defines the ini¬ 
tial value of S* at the moment of hard binary formation, 
expressed in the units of “full-loss-cone rate” ®, and 
the exponent v describes its decay with a. Theoretical 
models of the previous section and the results of Monte 
Carlo simulations suggest that in the collisionless case, 
/r = 0{1) and ly ~ 0.3 4- 0.6 in triaxial models (see Ta¬ 
ble [T]), while in the presence of relaxation ~ 0 and 
fj, I, scaling roughly as ^ ~ (A^*/10^)“^ in the spheri¬ 
cal and /i ~ (7V*/10®)“^/^ in the axisymmetric cases. 

As the binary hardens, GW emission becomes more 
and more effective. The instantaneous hardening rate 
due to GW is ^gw = l/(a7Gw)! where Tqw is given 
by Equation [51 We denote ogw to be the value of a at 
which 5'gw = 'S'*. Using the definitions of Oh m and 
'Sinfl ([5]) with A = 4, we obtain a simple estimate of the 
coalescence time for the case i/ = 0 (i.e. S'* = const): 


Ecoai 1.6 X 10*^ yr x 


X ^ 
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f'infi 

30 pc 

-1/5 
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( 22 ) 
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Let us now consider a qualitative model for the evolu¬ 
tion of the binary driven by both stellar-dynamical and 
GW hardening. 


d(l/a) 

dt 


= S*(a) -I- SGw(a, e) 


(23) 


= S 


*,h 


SGW,h 


/(e) 

a / /(eh) ’ 


where the values with subscript “h” refer to the moment 
of hard binary formation, and the eccentricity depen¬ 
dence /(e) is given by Equation [71 We further define a 
dimensionless parameter cgw = SGw,h/S*,h- In ah re¬ 
alistic situations, cgw ^ 1 ; meaning that at the early 
stage of evolution the hardening is driven by stellar en¬ 
counters and not by GW emission. If we assume that 
the eccentricity remains constant throughout the evolu¬ 
tion, then the above equation yields the following time 
to coalescence 


Ecoal 


1 r dy 
'S'*,hah y~’' + CGwy^ 

_^ 

4eGw 'S'*,hah V 5 + 1^ 5 -I- 


(24) 



Here 2 F 1 is the Gauss’ hypergeometric function and y = 
Oh/o. In the limit egw —t 0, the asymptotic value is 


Tcoai — 1-7 X 10® yr x 
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Fig. 7.— Evolution of the inverse semimajor axis (top panel) 
and eccentricity (bottom panel) for several triaxial models with 
no relaxation and the speed of light equal to 800 7V-body velocity 
units. Red dotted line is the Monte Carlo simulation, gre en d a shed 
line is the result of numerical integration of equations l|8 )l . (|10l> . 
and blue solid line is the reference Monte Carlo simulation without 
GW emission. Green dot marks the transition to GW-dominated 
regime (5* = S'gw)- The two evolutionary tracks ending at f 
1000 and t 1500 are for an equal-mass system {q = 1), and the 
other two are for mass ratio q = 1/9. 

The ratio between Oh and ogw; which describes how 
much the binary must shrink by stellar-dynamical pro¬ 
cesses before the GW emission takes over, is 

f \-TO 

OGw \30 pc/ \108M©y 

x[M/(e)]^ ((Y^)"^0.4U (26) 

The assumption of a constant eccentricity is not quite 
correct, though, as it both increases by stellar encoun¬ 
ters, as described by Equation [TUI and decreases due to 
GW emission, as given by Equation |8l In the final stage 
of evolution, when one may neglect 5** compared to S'gW) 
the following quantity is conserved: 

a-l (304 + 1216^)870/2299 (1 _ g2)-l ^ 

(27) 

Evolutionary tracks may be computed by numerical 
integration of the coupled system of equations (I8II0I23|) 
for a{t),e{t). Figure |7| shows several examples of evolu¬ 
tionary tracks computed using these equations, together 
with the results of Monte Carlo simulations of triaxial 
models with different initial eccentricity and mass ratio. 
We used the initial eccentricity as a free parameter in 
the evolutionary tracks, and adjusted to match the coa¬ 
lescence times from the Monte Carlo simulations: since 
the eccentricity varies rather erratically at early stages of 
evolution, it is hard to match these two curves without 
any tuning, but with this one free parameter one gets a 
quite good agreement at late stages of evolution. 

Figure [8] shows the coalescence times computed for 
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eccentricity 

Fig. 8 . — Coalescence time for triaxial models with Mbin = 
10® Mq and rjnfl = 30 pc, as a function of i nitial eccentricity. 
The initial hardening rate 5* is defined by CD with = 0.4 and 
ly = 1/3, in accordance with Monte Carlo simulations for an equal- 
mass binary. Green solid line is the result of numerical integration 
of the evolutionary track under our standard assumptions, purple 
dashed line is the evolutionary track computed witho ut s tellar- 
dynamical eccentricity growth (setting A = 0 in Eauation llOH . blue 
dot-dashed line is the track computed for a con stant 5* (i.e. z/ = 0), 
and red dotted line is the simple estimate lag which also assumes 
a constant hardening rate and neglects the changes in eccentricity. 

triaxial models with Mbin = 10® Mq, using the simple 
estimate (I22p and numerically computed evolutionary 
tracks, as a function of initial eccentricity. For e = 0, 
the coalescence time can be obtained analytically 
(Equation [25]) . For arbitrary eccentricity, it may be ap¬ 
proximated as 

T’coal « ITa? X (1 - e^) [fc + (1 - fc)(l - 6^)4] , (28) 

k = 0A + 0.1 logio(lWbin/10® Mq). 

It is instructive to compare the above expression to the 
simple estimate (1^ . which predicts Tcoai oc (1 — 
for the case of constant hardening rate and no evolu¬ 
tion of eccentricity. For e = 0 the latter underestimates 
the coalescence time by a factor ~ 3, but at higher ec¬ 
centricities the detailed evolutionary tracks come closer 
to the simple estimate, because the slowdown of stellar- 
dynamical hardening is compensated by the increase in 
eccentricity at the same evolutionary stage. The rather 
weak trend of the eccentricity-dependent factor with 
Mbin reflects the larger ratio between Oh and acw for 
smaller Mbin (Equation [551) , thus for them the stellar- 
dynamical increase in e is larger and consequently brings 
GW on stage earlier. 

The coalescence time very weakly depends on Mbin “ 
if we adopt the M, — a relation in the form (jH]), then 
Tcoai oc It also only moderately depends on the 

initial eccentricity: despite the much steeper dependence 
of the GW hardening rate on e, most of the time is spent 
on the stellar-dynamical hardening stage. As a conse¬ 
quence, the efficiency of stellar-dynamical hardening ^ 
is as important as the eccentricity. In the triaxial case, 
the efficiency is high enough and it decays slowly enough 
that for all reasonable parameters the coalescence time 
is shorter than the Hubble time. It also rather mildly 
depends on the fraction of chaotic orbits rj, which itself 
is determined by triaxiality; the caveat is that the latter 


changes in the course of evolution, but generally even a 
slight triaxiality is enough to support the required pop¬ 
ulation of centrophilic orbits. 

By contrast, in the collisionless axisymmetric case, the 
hardening rate slows down so rapidly that the binary 
stalls at a separation too large for efficient GW emis¬ 
sion, unless the eccentricity is very high. If we take into 
account relaxation-driven repopulation of the loss cone, 
then Equation|22|suggests that the coalescence time may 
be shorter than the Hubble time for Mbin 10® Mq and 
moderate eccentricity (under the optimistic assumption 
that the hardening rate, determined from Monte Carlo 

_ 2 ^ /2 

simulations to be proportional to W ' , stays at this 
relatively high value for a much longer time than these 
simulations were run). But in a realistic galaxy, even 
relatively minor perturbations from axisymmetry would 
create a sufficient reservoir of centrophilic orbits, whose 
draining maintains a much higher hardening rate than 
the relaxation can provide. 

5. EVOLUTION OF MERGER REMNANTS 

The isolated models considered above serve as a con¬ 
trolled experiment that helps to understand the physical 
mechanisms responsible for the joint evolution of stars 
and the binary. However, the initial conditions were 
somewhat artificial in that we have set up initial models 
in almost perfect equilibrium. In the cosmological set¬ 
ting, binary SBHs are expected to form via galaxy merg¬ 
ers, and the merger remnants could well have a complex 
and evolving structure quite unlike our idealized models. 

In this section we consider a limited se t of me rger sim¬ 
ulations, similar to those of iKhan et al.l (120111) . We set 
up two identical spherical 7 = 1 Dehnen models with 
scale radius and mass equal to unity, each containing a 
central SBH with mass M, = 0.01. They are put on an 
elliptical orbit with initial separation 20 and relative ve¬ 
locity 0.1; the first encounter between the galaxies occurs 
at t ~ 80, the second at t ~ 100, and by t = 110 the two 
nuclei are well merged and the binary is formed. The 
spherically averaged density profile of the central region 
of the merger remnant right after hard binary formation 
is well described by a 7 = 0.5 Dehnen profile with unit 
scale radius and total mass ~ 1.7. We evolve the system 
until t = 300, using the same direct-summation code 
((>GRAPEch, and extract a snapshot at a time t ~ 120 
when the binary semimajor axis a = 0.01 < ah', this 
snapshot then constitutes the initial conditions for the 
Monte Garlo simulations. We ran four simulations with 
particle numbers N = {128,256,512,1024} x 10®. The 
morphology of the merger remnants was roughly the 
same, but the initial eccentricity of the binary at the 
moment of its formation was systematically larger for 
higher-resolution simulations (FigurejO]), even though the 
orbit parameters of the merging galaxies were identi¬ 
cal; thus these four Monte Carlo models are not exactly 
equivalent. 

Unlike the isolated models, merger remnants do not 
exhibit genuine triaxial symmetry, but only a reflection 
symmetry (i.e. even after a rotation that aligns the ma¬ 
jor axis with the a;-coordinate axis, they are not in¬ 
variant with respect to a flip about either the x or y 
axes, but stay the same under a simultaneous inversion 
of both axes - like a barred spiral galaxy in which the 
spiral arms break the triaxial symmetry of the bar). We 
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Fig. 9. — Binary ecccentricity in merger simulations: solid lines 
are A'^-body and dashed lines are Monte Carlo models with relax¬ 
ation; from bottom to top: N = {128,256,512, 1024} X 10®. 

therefore kept all reflection-symmetric (even I and all m) 
terms in the spherical-harmonic expansion of the poten¬ 
tial in the Monte Carlo code. The density profile rotates 
with a non-uniform angular speed, which precludes the 
use of a rotating reference frame; since the rotation is 
rather slow, we just updated the expansion coefficients 
frequently enough (each 1 time unit) to track this figure 
rotation. We only used terms up to quadrupole (1 = 2), 
therefore smoothing out all irregularities and small-scale 
structure that might be present in the merger remnant, 
and retaining only the global non-spherical features. We 
have checked that using higher-order harmonics does not 
substantially change the results. In addition, we have 
run simulations with imposed axisymmetry (setting all 
TO ^ 0 terms to zero). For each of the four values of N, 
we have performed Monte Carlo simulations with relax¬ 
ation (corresponding to AC = A^ in the actual A^-body 
system) and without relaxation. 

The results, shown in Figure 1101 can be summarized 
as follows. In the A^-body simulations, the hardening 
rate was found to be almost independent of N, in agree¬ 
ment with other studies. However, the highest-A^ model 
had a slightly lower hardening rate at late times. The 
weaker A'^-dependence of hardening rate in merger rem¬ 
nants compared to idealized isolated models could have 
a number of reasons, e.g. perturbations from the decay¬ 
ing large-scale clumps and inhomogeneities in the merger 
remnant are presumably independent of N and add up to 
the conventional two-body relaxation. This conjecture is 
supported by the fact that the rate of energy and angu¬ 
lar momentum diffusion measured in the simulations is 
almost independent of N and much higher than the rate 
expected from two-body relaxation. The eccentricity, be¬ 
ing quite high at the formation time, increased further 
by the end of the simulation, up to > 0.98 for N = 10®. 
We note that the measured hardening rate was ^ 3 times 
lower than 5'fun, again underlining that even in strongly 
asymmetric and dynamic systems the loss cone is at least 
partially depleted at the highest binding energies. 

While the results Monte Carlo simulations in this sec¬ 
tion should be regarded as preliminary, it is remarkable 
that the agreement in hardening rates between A^-body 
and Monte Carlo simulations that included relaxation 



0 50 100 150 200 250 300 


time 



time 


Fig. 10.— Evolution of inverse semimajor axis of the binary in 
merger simulations. Top panel compares A^-body (solid lines) and 
Monte Carlo models with relaxation (dashed lines), from left to 
right: N = {128, 256, 512,1024} X 10® (curves are shifted horizon¬ 
tally for clarity). Bottom panel shows the long-term evolution of 
collisionless Monte Carlo models, with imposed axisymmetry (dot¬ 
ted lines) and with only reflection symmetry (i.e., nearly triaxial 
models with figure rotation, dashed lines). Solid line shows the 
average hardening rate in collisional simulations. 

is reasonably good (Figure (TUI top panel). Relaxation 
driven by large-scale fluctuations in the merger remnant 
might be underestimated in the Monte Carlo method 
because of two factors: (a) the potential represented 
by spherical-harmonic expansion smooths out small-scale 
clumps, (b) the interval between potential updates is 
short enough to track global changes in the potential, but 
may not represent higher-frequency transient perturba¬ 
tions. Thus we should expect somewhat lower hardening 
rate in Monte Carlo simulations compared to V-body 
simulations; in fact it turned out to be comparable and 
even sometimes higher. We defer a more detailed study 
of Monte Carlo models of merger remnants for a future 
work. 

In the collisionless regime, the hardening rate is ini¬ 
tially only slightly lower than in the V = 10® simulation 
with relaxation. However, it slows down later on, much 
like in the case of isolated triaxial systems (Table [IJ. 
Again, this is due to two factors: depletion of centrophilic 
orbits and gradual decrease of triaxiality. The shape of 
the merger remant is moderately flattened (z/a; ~ 0.75 
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at all radii throughout the simulation), and the triaxi- 
ality is quite subtle: yjx ~ 0.9 right after merger and 
increases to > 0.97 toward the end of simulation. Never¬ 
theless, this appears to be enough to sustain the reservoir 
of centrophilic orbits needed to keep the binary shrink¬ 
ing. Even with zero eccentricity, the binary would merge 
in ^ 0.5 Gyr; with such high eccentricity as in our simula¬ 
tions, this time would be an order of magnitude shorter. 
If we force the potential to be axisymmetric, the evolu¬ 
tion slows down much faster iFigurelTOl bottom panel); 
by the end of simulation a ~ O.Oboh, which is not enough 
to ensure merger in a Hubble time unless the eccentricity 
is > 0.8. 

6. DISCUSSION AND CONCLUSIONS 
6 .1. Summary of our results 

In the present work, we have considered a stellar- 
dynamical solution to the final-parsec problem - the evo¬ 
lution of binary SBHs driven by encounters with stars in 
a galactic nucleus, as its orbit shrinks from the radius of a 
hard binary (oh 1 pc) to the separation ogw ~ 10 “^ Oh 
at which GW emission becomes effective. The central 
difficulty is to ensure a continuous supply of stars onto 
low-angular-momentum orbits, where they can be scat¬ 
tered by the binary and carry away its energy and an¬ 
gular momentum. This region of phase space, dubbed 
the loss cone, is quickly depleted by the binary once 
it becomes hard, so an efficient refilling mechanism is 
required to enable continued hardening. In a perfectly 
spherical galaxy, this can only be achieved by two-body 
relaxation, which is not sufficient to bring the binary 
to coales cence in a Hubble t ime, except in the smallest 
galaxies (jMerritt et al.l[200^ - hence the problem. We 
have focused on the additional mechanisms of loss-cone 
repopulation that exist in non-spherical (axisymmetric 
and triaxial) galaxies. 

We addressed this problem using a variety of meth¬ 
ods, but primarily with Raga 13 a novel stellar-dynamical 
Monte Carlo code that is able to follow the evolution of 
non-spherical stellar systems under the influence of two- 
body relaxation, the magnitude of which can be adjusted 
- unlike conventional iV-body simulations in which the 
relaxation is essentially determined by the number of par¬ 
ticles. We extended the code to include interactions be¬ 
tween stars and the massive binary, by following the tra¬ 
jectories of particles in the superposition of the smooth 
galactic potential plus the time-dependent potential of 
the two SBHs as they orbited one another, and used con¬ 
servation laws to determine the reciprocal changes in the 
binary orbital parameters. 

We used a large suite of direct-summation 7V-body sim¬ 
ulations to verify that the Monte Carlo code accurately 
describes the evolution of the binary for a wide range of 
parameters - mass ratio, eccentricity, number of parti¬ 
cles and the shape of the galaxy model. Then we ex¬ 
tended the Monte Carlo simulations into the range of fV* 
much larger than is presently accessible for conventional 
A^-body simulations, including the intriguing collisionless 
limit (Ni, = oo), which is nearly achieved in real galaxies. 
We determined the scaling laws and asymptotic behav¬ 
ior of the stellar-dynamical hardening rate from rather 

^ The code is available at http://td.lpi.ru/-eugvas/raga/ 


simple analytic arguments, and confirmed these findings 
with the Monte Carlo simulations. Taking into account 
GW emission, the evolution of binary semimajor axis a 
and eccentricity e can be described by a simple system of 
differential equations, which we used to determine the co¬ 
alescence time as a function of initial parameters of the 
binary and the galaxy; these evolutionary tracks were 
again verified by Monte Carlo simulations for a few test 
cases. 

Our main results can be summarized as follows: 

1. The binary continues to shrink as long as there are 
stars in the so-called loss region - the region of 
phase space from which stars can precess into the 
loss cone proper (that is, < Llc = 2 GMbina) 
due to non-spherical torques. In the spherical case, 
the loss cone and loss region are the same. In 
the triaxial case, the loss region consists of mostly 
chaotic orbits, whose total mass is roughly propor¬ 
tional to the mass of galaxy with a coefficient r; <C 1 
determined by the degree of flattening and triaxi- 
ality, and unless the galaxy is nearly axisymmetric, 
the mass of stars in the loss region is much larger 
than Mbin. In the axisymmetric case, the volume 
of the loss region shrinks with a as y/ija, halfway 
between the spherical and triaxial cases. 

2. The stellar-dynamical hardening rate S'* is always 
smaller than the value Sfuii corresponding to a full 
loss cone. This is explained by a gradual depletion 
(draining) of the loss region, which occurs faster 
at high binding energies. As a rough estimate, de¬ 
creasing a by a factor of two requires ejection of 
stars with total mass ^ Mbin- In the collisionless 
limit, the hardening rate declines with a very slowly 
for a triaxial system, because the total mass of loss 
region stars is typically large compared to Mbin- 
By contrast, in an axisymmetric system, the vol¬ 
ume of the loss region also shrinks with a, and it 
depletes much faster; thus the hardening rate de¬ 
clines rapidly and drops nearly to zero at a ;< O.loh- 
This is not enough to bridge the gap to the GW- 
dominated regime. 

3. Taking into account relaxation-driven repopulation 
of the loss region does not change the results in the 
triaxial case very much, because the hardening is 
dominated by draining of stars that are initially 
in the loss region, but accounting for relaxation 
does change the dynamics in the axisymmetric and 
spherical cases dramatically. After the initial pop¬ 
ulation of the loss region is nearly depleted, it is re¬ 
filled by two-body relaxation, maintaing the hard¬ 
ening rate at a nearly constant level that depends 
on TV*. The axisymmetric case offers a didactic 
example of how much a system with W ^ 10®, 
typical of present-day, high-fidelity Wbody simu¬ 
lations, or even TV* ^10®, can differ from the col¬ 
lisionless limit in its long-term evolution. 

4. Coalescence times estimated for triaxial galaxies 
weakly depend on the mass of the binary, its mass 
ratio, or the degree of triaxiality (provided that the 
departure from axisymmetry is larger than a few 
percent). Coalescence times fall in the range from 
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a few Gyr for almost circular binaries, to < 10® yr 
for very eccentric ones. For a given combination 
of the binary mass Mbin and its radius of influence 
Tinfl, the coalescence time can be computed using 
Equation [52] for a circular orbit, and using Equa¬ 
tion |5S| for an eccentric orbit, where the typical 
values of the dimensionless parameters /r ^ 0.2 1 

and V ^ 0^0.5 are listed in TablejTJ This time is up 
to a factor of few times greater than the simple es¬ 
timate (1551) that does not account for the decrease 
of the hardening rate with time. 


6 .2. Comparison with previous work 

Consider first the spherical case, which has been 
extensively studied by IV -b ody simulations (e.g. 
Quinlan fc Hernauisd 119971: SElosavlievic fc Merritti 


20011: Hemsendorf et ahl 12002 : iBerczik et al.l 120051 : 

Merritt et al.l 200711. methods based on scattering exper- 


iments ( Quinlanlll996l : lSesanal201(ilMeiron fc Laoill26l2L 

e.g.), or Fokker-Plan c k mo dels (|Milosavlievic fc MerrittI 
I2003bl : IMerritt et al.1 120071) . It is generally accepted 
that in the spherical case the binary quickly depletes 
the loss cone and its evolution nearly stalls at a value 
of a just a few times smaller than Oh. A number of 
effects may moderately decrease the stalling radius or 
increase the relaxation rate. The secondary slingshot - 
the re-ejection of stars that have once interacted with 
the binary but did not gain enough energy to escape 
the galaxy - lead s to a gradual (oc logt) incre ase of 
1 /a at late times ([Milosavlievic fc Merritt! l2003bD . and 
this trend was indeed found in our collisionless Monte 
Carlo simulations. Once these stars are completely 
eliminated, evolution o f the binary fi nally stalls in the 
pure^ collisionless case: IMerrittl (I2006D and lSesana et al.l 
(120071) found that the stalling radius is only a few times 
smaller than Oh over a wide range of binary mass ratios 
and cusp density profiles. Similarly, time-dependent 
solution of the Fokker-Planck equation describing the 
relaxation of stars in angular momentum yields a higher 
rate of loss-cone repopulation at early times than the 
steady -state flux, due to sharper grad ients in the phase 
space ([Milosavlievic fc Merrittll2003bD : this is taken into 
account automatically in the Monte Carlo scheme, and 
does not substantially affect the overall evolution. 

Brownian motion is not accounted for in our method, 
and this may be the reason for the discrepancy between 
Monte Carlo and A^- body hardening rates at hig h N 
in spherical systems. [Quinlan fc Hernauisd (119971) and 
iChatteriee et al.l ([2003D argued that wandering of the bi¬ 
nary may explain the very weak dependence of harden¬ 
ing rate on N found in their simulation s. However, as 
summarized bv IMakino fc Funatol ([2004D . several factors 
complicate the interpretation of their result: the A^-body 
algorithm used in that paper is less collisional than tra¬ 
ditional A^-body schemes, but not entirely free of relax¬ 
ation, and the unequal masses of particles mean that the 
granularity of the potential depends on the evolution¬ 
ary stage. The amplitude of Brownian motion of the 
binary, while larger than for a single SBH of the same 
mass ([MerrittI 120011) . scales roughly as N and in 
real galaxies would be smaller than the size of the loss 
cone for most part of the evolution. Using somew hat 
different arguments, iMilosavlievic fc MerrittI (l2003bD es¬ 
timated the timescale of loss-cone refilling by Brownian 


motion and concluded that this effect is unlikely to sub¬ 
stantially affect the evolution. 

Consider next the case of non-spherical galaxies, which 
seem to offer a more promising way to solv e th e final- 
parsec problem via collisionless dynamics. [^ (120021 ) 
estimated draining rates for the loss regions in axisym- 
metric and triaxial galaxi es, using a rguments similar to 
those in our section l43l IYuI (I2002D concluded that the 
loss region in triaxial galaxies is not likely to be depleted 
if the flattening parameter, responsible for the fraction of 
centrophilic orbits, is > 0.05, although the plots in which 
she showed evolution timescales demonstrate the effect of 
gradual depletion of the loss region only for the c ase o f 
small flattening. For axisymmetric systems, lYul ([2002D 
estimated that the loss region ('or “loss wed ge”, in the 
terminology of iMagorrian fc Tremainelll999l) can be de¬ 
pleted rather quickly in many cases, and the relaxation- 
limited evolution timescales are still longer than the Hub¬ 
ble time. Thus our conclusions qualitatively agree with 
t hat study. _ 

IMerritt fc PoonI (I2004D considered self-consistent tri¬ 
axial models of galactic nuclei with SBHs which had a 
significant fraction of centrophilic orbits. They solved the 
evolutionary equations similar to ([121151) , and found that 
the hardening rate is nearly independent of time, but 
scales with the square of the fraction of chaotic orbits. 
Analysis of the data plotted in their figure 6 suggests 
that 1/a oc ^0-7—0-9^ more in li ne with our findings in the 
present study and injPapeO- For the case of a singular 
isothermal cusp (7 = 2 ), Equation (ITTl) suggests that the 
asymptotic hardening rate scales as the squared fraction 
of chaotic orbits rj, in agreement with their results. Their 
equation 55 implies a hardening rate ~ jy^iS'fun, and our 
Monte Carlo simulation of a 7 = 2 Dehnen model has 
on average a similar hardening rate with rj = 0.2, even 
though it declines with time. Thus the basic conclusion 
of that paper, that even a moderate amount of triaxial- 
ity is sufficient to drive the binary to coalescence in less 
than a Hubble time, is corroborated by our simulations 
and asymptotic analysis, even though some details differ 
(most importantly, their neglect of the slowdown of the 
hardening rate). 

Most other studies to date have assumed or in¬ 
ferred that the loss cone must be kept nearly full 
in non-spherical systems, using vario us arguments. 
iHollev-Bockelmann fc SigurdssonI ([2006D explored the 
properties of orbits in a triaxial galaxy with a central 
SBH. They estimated that the time required to change 
the angular momenta of particles due to collisionless 
torques was much shorter than the Hubble time, and 
concluded that the loss cone must remain full. However, 
this argument does not take into account the draining of 
the loss region, and they did not analyze the evo lution of 
their m odel under this process. More recently, ILi et al.1 
([2014D performed a similar analysis for a nearly axisym¬ 
metric model, which in fact was slightly triaxial in the 
central part. They computed the mass of stars belonging 
to orbits that are able to come into the loss cone with 
a size corresponding to the initial stage of binary evo¬ 
lution, which was several times larger than Mbin- From 
this they concluded that the loss cone should remain well 
populated during the subsequent evolution, but their es¬ 
timate did not take into account that the volume of the 
loss region also shrinks along with the binary semima- 
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jor axis in the axisymmetric case. It is unclear whether 
the slight triaxiality of their model would be sufficient 
to support enough truly centrophilic orbits during the 
e ntire evo l ution. 

iSesanal (|201flt l considered a hybrid model for bi¬ 
nary evolution, based on the hardening and eccentricity 
growth rates computed from scattering experiments. He 
assumed that after the initial phase of formation of a 
hard binary, accompanied by the erosion of the stellar 
cusp at r < rinfl, the subsequent evolution occurs in the 
full-loss-cone regime, i.e. the hardening rate is given by 
‘S'uniform (Equation [2|) with density and velocity disper¬ 
sion computed at rinfl. Thus his evolutionary tracks are 
similar to our calculations in Section 14.41 with parame¬ 
ters /i = 1 (the loss cone is full) and ly = 0 (the harden¬ 
ing rate does not decrease with time). As Figure [ 8 ] and 
Equation (j22p show, in this case the coalescence time is 
shortened by a factor of few with respect to the more 
conservative assu mptions propose d in ou r study. 

More recently, iSesana fc KhanI (1201511 compa red the 
evolut ionary tracks from the hybrid model of ISesanal 
(|2010l l wi th those obtained by A^-body simulations of 
mergers (jKhan et al.l 120121 1. They found reasonable 
agreement for the eccentricity evolution and hardening 
rate, even though the latter was somewhat lower and de¬ 
clined with time in A^-body simulations. They ascribed 
this to the gradual decrease of the density and corre- 
spo nding increase of rinfl. As we have argued in Sec¬ 
tion 14.31 the decline of hardening rate is mostly caused 
by the depletion of centrophilic orbits at all radii, rather 
than simply the decrease of density at the influence ra¬ 
dius. For instance, in our collisionless simulations of 
isolated models with 7 = 1 ( 2 ), the actual hardening 
rate dropped by a factor of 5 (lO) by the end of simula¬ 
tions, while the density at rinfl, and correspondingly the 
full-loss-cone hardening rate Sfuii, decre a.sed by less than 
30%. The coalescence times quoted in ISesana fc KhanI 
(|2015ll are longer than ours due to a different Mbin — rinfl 
relation adopted in that paper. We stress, however, that 
their estimates are based on the hardening rates from col- 
lisional A^-body simulations, while coalescence times for 
collisionless systems, advocated in the present study, are 
up to a few times longer for the same galaxy parameters. 

Studies based on A-body simulations have gener¬ 
ally observed little or no dependence of the hardening 
rate on A in non- spherical galaxy models that were 
forme d via mergers (|Preto et al.ll2011 |: |Khan et al. 20H, 
2 OI 2 II or created as isolated models (jBerczik et al. I 2 OO 6 : 
Khan et al.ll2013ll . This result has been interpreted as an 
indication tha t the lo ss cone remains nearly full, although 
iBerczik et al.l (|2006fl reported that the hardening rate 
was gradually decreasing with time, suggesting that the 
reservoir of centrophilic orbits was being depopulated. 
A si milar trend can b e seen in other merger simu lations 
(e.g. lPreto et al.l[201lL Fig.l, or lKhan et al.ll201^ Fig.2); 
in the latter paper, the models with steeper cusps dis¬ 
played a systematically more rapid decline of hardening 
rate with time. All these trends are in agreement with 
our findings, even though the previous studies did not 
highlight them. We note, however, that in our simula¬ 
tions the hardening rate always turns out to be substan¬ 
tially lower than S'fun for large enough A. This might 
seem to be in contrast with other studies that report 
a hardening rate comparable to S'fun for non-spherical 


systems. However, a more detailed examination suggests 
that the apparent discrepancy can be attributed to differ¬ 
ent definitions of the full-loss-cone rate, and to different 
no rmalizations of the stellar density profile. For instance, 
in iKhan et al.l (j2013[l and iHollev-Bockelmann fc KhanI 
(j201,5h the flattened models were created by adiabatic 
squeezing of the original density profile, and hence their 
scale radii are roughly a factor of two smaller than ours, 
as can be seen in Figure 1 of the latter paper. From Equa¬ 
tion ([5]), it is apparent that this translates to a hardening 
rate ^ 6 times higher than ours. We have re-simulated 
their models with our A-body code and found generally 
good agreement with their results. 

Most importantly, as we have argued in iPaoer II and 
this study, the hardening rates in A-body simulations 
are dominated by, or at least have a significant contribu¬ 
tion from, collisional effects even for A ~ 10®, thus it is 
not easy to extrapolate these results to real galaxies. Us¬ 
ing the Monte Carlo method, we were able to reach the 
collisionless regime, which turned out to be very different 
for axisymmet ric and triaxial galaxies, while in A-body 
simulations of iPaoer^ they looked nearly the same. 

6.3. Single and binary black holes 

It is also instructive to compare loss-cone theory in the 
single and binary SBH cases. The first obvious difference 
is the much larger size of the loss cone in the case of a 
binary. As a consequence, the relaxation-driven repopu¬ 
lation of the loss cone almost always occurs in the empty- 
loss-cone regime; on the other hand, collisionless changes 
in angular momentum due to non-spherical torques occur 
on the same dynamical timescale as the depletion of the 
loss cone, so that it always remains partially populated in 
non-spherical systems. A second important factor is that 
the size of the loss cone decreases as the binary shrinks, 
while in the case of a single SBH it can only grow. Third, 
the mass of stars needed to be delivered into the loss cone 
of the binary is a few times larger than the mass of (the 
lighter component of) the binary, while for single SBHs 
the accreted mass in stars is typically small compared to 
the SBH mass. 

These three factors explain the fundamental differ¬ 
ence between collisionless spherical and axisymmetric 
systems, on the one hand, and triaxial ones, on the other 
hand. For the latter ones, the evolution is almost en¬ 
tirely driven by draining of centrophilic orbits, whose 
total mass is much larger than M, and furthermore does 
not depend on the size of the loss cone. Interestingly, in 
the case of a single SBH most of captured stars arrive 
from regular pyramid orbits inside rinfl, while in the case 
of the binary the loss region consists mainly of chaotic or¬ 
bits outside rinfl; this explains the slightly different time 
dependence of the draining rates. In the axisymmetric 
case, the volume of the loss region composed of chaotic 
orbits that can be delivered into the loss cone by colli¬ 
sionless torques, shrinks along with the binary semimajor 
axis, and its orbit population is nearly depleted before 
the binary reaches the GW-dominated regime. The sub¬ 
sequent evolution is determined by the rate at which this 
loss region is repopulated by relaxation. Since the vol¬ 
ume of this region is still much larger than the volume of 
the loss cone proper, it is more easily repopulated in ax¬ 
isymmetric than in spherical systems. The same is true 
for single SBHs; the fact that for them the difference 
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between axisymmetric and triaxial systems is much less 
than b etween spherical and axisymmetric ones (|Vasilievl 
l2014al Figure 4) stems largely from the adopted isotropic 
(non-depleted) initial conditions for the relaxation. On 
the other hand, in the case of a massive binary the phase 
space is already depleted out to much larger values of 
angular momentum than the current loss cone bound¬ 
ary, thus it takes longer for the relaxation to resupply 
the loss region. 

In short, loss cone theory in non-spherical systems is 
a delicate interplay between collisional and collisionless 
effects, and the outcome depends on the evolutionary his¬ 
tory of the loss cone, as well as the changes in the global 
structure of the system (its shape and phase-space gradi¬ 
ents) . Only using a combination of various approaches - 
7V-body simulations, orbit analysis, Monte Carlo meth¬ 
ods and scaling arguments - can one hope to understand 
the behavior of realistic stellar systems. 


6.4. Conclusions 

The evolution of binary SBHs in gas-poor galaxies is 
determined by the rate of slingshot interactions with 
stars in the loss cone - the low-angular-momentum re¬ 


gion of the phase space. The fact that the loss cone is 
quickly depleted in idealized spherical systems gave rise 
to the final-parsec problem. Repopulation of the loss 
cone occurs both due to collisional and collisionless ef¬ 
fects; the latter are only relevant in non-spherical sys¬ 
tems. We have developed a Monte Carlo method that 
can efficiently deal with both collisionless and collisional 
evolution, and used it to show that in the collisionless 
limit, the repopulation is efficient if the galaxy is even 
slightly triaxial. To the extent that mergers result in 
galactic shapes that are not exactly axisymmetric, our 
results imply that the final-parsec problem does not ex¬ 
ist in most galaxies. 
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fruitful discussions. We are grateful to the referee for 
valuable comments. This work was supported by the 
National Aeronautics and Space Administration under 
grant no. NNX13AG92G to D.M. We acknowledge the 
use of computing resources at GIERA funded by NSF 
PHY-1126812. 
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